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ABSTRACT 

Interval  integration  is  used  to  obtain  inclusions  of  integral  operators  of 
the  form 

g(u)  ( 

which  can  be  carried  out  on  a  computer.  The  resulting  inclusions,  combined  with 
interval  iteration,  are  used  to  compute  guaranteed  upper  and  lower  bounds  for 
solutions  of  integral  equations  of  the  form 

u  -  g(u) 

j  ,  ‘  -  Vl’  /  '1 

for  86s.  It  is  also  possible  to  establish  existence  or  nonexistence  of 
solutions  of  integral  equations  in  given  regions  on  the  basis  of  results  of  the 
computation.  Examples  of  applications  of  this  technique  to  linear  and  nonlinear 
integral  equations  are  eigenvalue  problems  for  linear  integral  operators  are 
given. 
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Key  words:  Integral  equations.  Eigenvalue  problems.  Error  bounds.  Interval 
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SIGNIFICANCE  AND  EXPLANATION 


This  report  is  the  text  of  an  invited  address  to  be  given  at  the  11th  IMACS 
World  Congress,  to  be  held  in  Oslo,  Norway,  on  August  5-9,  1985.  The  theories 
of  interval  integration,  interval  iteration  and  their  applications  to  integral 
equations  have  been  developed  in  previous  papers.  This  report  surveys  these 
results,  but  concentrates  on  computer  implementation  by  discussing  various  ways 
in  which  inclusions  of  integral  transforms  can  actually  be  represented  on  a 
computer.  The  computations  based  on  these  inclusions  are  self -validating,  in 
that  existence  of  solutions  and  upper  and  lower  bounds  for  the  solutions  of 
linear  and  nonlinear  integral  equations  and  eigenvalue  problems  caui  be 
guaranteed  on  the  basis  of  computed  results  alone.  It  is  assumed,  of  course, 
that  the  computing  environment  supports  interval  computation.  Since  support  for 
interval  arithmetic  is  now  required  by  the  IEEE  standard  for  floating-point 
arithmetic,  and  provided  by  Pascal-SC  for  microcomputers  and  ACRITH  for  IBM  370 
mainframe  computers,  the  need  to  make  special  provision  for  interval  arithmetic 
is  only  necessaury  on  outdated  systems. 

Some  simple  examples  illustrating  the  techniques  presented  in  the  paper  aure 
cited. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 
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COMPUTABLE  BOUNDS  FOR  SOLUTIONS  OF  INTEGRAL  EQUATIONS 


L.  B.  Rail 


Summary 

Interval  Integration  is  used  to  obtain  inclusions  of  integral  operators  of  the  form 

g(u) (s)  -  JTg(s,t,u(s) ,u(t) )dt  (1) 

which  can  be  carried  out  on  a  computer.  The  resulting  inclusions,  combined  with  interval 

iteration,  are  used  to  compute  guaranteed  upper  and  lower  bounds  for  solutions  of  integral 

* 

equations  of  the  form 


u  -  g<u)  (2) 

for  s  €  S.  It  is  also  possible  to  establish  existence  or  nonexistence  of  solutions  of 
integral  equations  in  given  regions  on  the  basis  of  results  of  the  computation.  Examples 
of  applications  of  this  technique  to  linear  and  nonlinear  integral  equations  and  eigenvalue 
problems  for  linear  operators  are  given. 

Integral  Equations 

The  integral  equation  expressed  by  (2)  includes  many  of  the  varieties  of  integral 
equations  of  importance  in  applications.  For  example,  for 
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■v"  .. 


g(s,t/U(  ’.)  ,u(t)  )  -  f<X)  +  AK(s,t)u(t) , 


the  result  Is  the  linear  Integral  equation  of  Fredholm  type  and  second  kind. 


u(s)  -  f(x)  +  x/TK(s,t)u(t)dt, 


while  the  choice 


g(s,t,u(s) ,u(t) )  -  f (s,t,u(t) ) , 


gives  the  nonlinear  integral  equation  of  Urysohn  type 


u(e)  -  J  f(s,t,u(t))dt. 


and  so  on  [12] .  In  general,  S  and  T  are  subsets  of  some  possibly  high-dimensional  real  or 
cosgilex  spaces.  For  definiteness,  we  will  take  S  *  T  -  [0,1],  and  consider  integral 
equations  of  the  form 


u(s)  -  J|lg(8,t,u(8)  ,u(t)  )dt, 


since  the  discussion  given  here  will  carry  over  immediately  to  the  more  general  case.  It 
is  also  worth  noting  that  (7)  includes  Volterra  integral  equations  as  the  special  case 
g(s,t, v,w)  =  0  for  t  >  s. 


Interval  Analysis 


Our  purpose  is  to  find  functions  u(t),  u(t)  such  that 


u(s)  <  u(s)  <  u(s) 


.*•  .**  .' 

.  • .  »  .  ■  .*  .  • . 


v,.' 

'■.*  ■.*  v  /v 


v*  V  *,•  V  •  J 

.V-V-V-  Si 

•'  .  **.  *  «  •* 


0  <  a  <  1,  in  other  words,  an  interval-valued  function 


0(s)  -  [u ( e ) ,  u(e)] 


such  that  u(s)  €  U(s)  for  s  €  [0,1]  [1]-  The  process  described  in  [12]  for  this  purpose  is 
a  combination  of  interval  integration  [1]  and  interval  iteration  [9]. 


Interval  Integration 


The  integral  of  an  interval  function  U(t)  is  defined  to  be 


/Jo(s)ds  -  [/’u(s)ds,  /’u(s)ds]. 


where  under  and  overbars  on  the  integral  signs  within  the  brackets  denote  respectively 
lower  and  upper  Darboux  integrals  [1].  Since  these  Darboux  integrals  always  exist,  it 
follows  that  all  interval  (and  hence  all  real)  functions  are  integrable  in  this  sense. 


Interval  Iteration 

Interval  iteration  far  the  fixed-point  problem  (2)  works  like  this  [9],  [12]:  Suppose 
that  one  has  an  Interval  inclusion  G  of  the  operator  g,  that  is,  an  operator  from  interval 
functions  to  interval  functions  such  that  g(u)  t  G(U)  for  each  u  t  U.  Then,  starting  with 
an  intial  interval  function  UQ,  one  computes 


Vi  “  °n  nG<V' 


Interval  iteration  has  the  following  properties  [9] : 

1.  If  u  €  Uq  is  a  solution  of  (2),  then  the  interval  iteration  will  converge  to 


0*  -  n.  U  , 

n«D  n 


*.  ",  \  *,  •,  », 


" •  ■*  ,*  ,•  *#e  *_»  .  • 


*  •.*  a *  .  V  ",  •- 


and  u  €  0* 


2.  If  for  some  integer  n. 


Un  n  G(0„)  -  0.  (13) 

the  empty  set,  then  there  are  no  solutions  u  of  (2)  in  the  initial  interval  UQ. 

3.  If 

G(U„)  C  U„  (14) 

for  some  integer  n,  then  a  solution  u  of  (2)  exists  in  the  initial  interval  UQ. 

Thus,  interval  iteration  can  be  used  to  prove  existence  of  solutions  if  (14)  holds, 
nonexistence  by  (13),  or  to  obtain  improved  bounds  for  solutions  by  ( 12) .  The  problem  in 
actual  computation  is  then  to  construct  an  interval  inclusion  G  of  the  integral  operator  g 
in  (1).  One  way  to  do  this  is  to  construct  an  inclusion  G(s,t,U(s),U(t) )  of  the  integrand 
by  using  interval  arithmetic  [6] ,  and  then  use  interval  Integration  to  include  the 
operator.  Here,  by  "interval  arithmetic",  we  mean  interval  inclusions  of  the  ordinary 
library  functions  available  on  a  computer  in  addition  to  the  arithmetic  operations 
themselves.  He  now  present  several  ways  to  implement  such  a  construction  on  an  actual 
computer.  It  is  assumed,  of  course,  that  interval  computation  is  properly  supported,  as  in 
Pascal-SC  or  ACRITH. 

Inclusion  of  Integral  Operators 


Interval  Arithmetic 

A  straightforward  method  to  construct  an  inclusion  of  the  integral  operator  ( 1 )  simply 
uses  interval  arithmetic.  Here,  the  interval  S  ■  T  -  (0,1]  is  partitioned  into  n 
subintervals 


-4- 


St  m  T^  * 


with  endpoints  s^.,  <  s^,  sQ  *  °»  an  *  1#  and  widths  wlS^)  *  Sj-s^.j.  Further,  we  define 
an  interval  step  function  U  by  U(e)  «  0^«  J  for  s  €  S^«  Then,  if  u  €  0,  it  follows 

that 

g(s,t,u(s),u(t))  t  G(Si,Tj,Oi,U3)  (15) 

for  s  €  S^,  t  €  T y  Thus, 

n 

U(s>  i  l  G  (0  ,0  )*w(T  ),  (16) 

j.1  13  *  3  3 

where  G1j(Ui,U^)«G(St,T^,U1,U^),  because  the  Interval  Rlemann  sum  on  the  right  Bide  of  (16) 
contains  the  interval  integral  of  G(a,t,U(s) ,U(t) >  11].  In  this  case,  the  interval 
iteration  (11)  can  be  carried  out  in  the  form 

°T  ■  0>ij(«*’-(V  + 

(17) 

n 

♦  l  Glj<°J,®^,*w(Tj>]  n 

a  "Gauss-Seidel"  iteration  which  makes  use  of  updated  Information  as  it  is  produced,  and 
has  the  sans  convergence  properties  as  (11)  (13].  This  simple  approach  to  solving  integral 
equations  has  the  following  advantages i 

1.  It  is  easy  to  implement  on  a  computer. 

2.  No  special  properties  of  g(s,t,v,w)  such  as  positivity  or  monotonicity  are 

-5- 


.'-V* 


.*■  '.'A 'V-V-V-V*’ 

■  ’»  I  VI*'  I  I  <  •  .  •'  »  -*  »  -  ■  •  a 

••  •  ,*  ..  •  •  l.  *  .*»  V*  . /  ./ 


required. 
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3.  The  computation  is  self-validating,  with  existence  of  a  solution  of  the  integral 
equation  and  lower  and  upper  bounds  for  it  guaranteed  if  (14)  holds. 

4.  A  form  of  "superconvergence"  holds  at  the  nodes  8i  for  i  -  1,...,n-1.  since  u(8i) 
{  Djf  and  u(s^)  (  o"+1,  it  follows  that 

u‘8i>  t  uj  n  oJ+1,  (18) 

which  results  in  narrower  bounds  for  u ( s )  at  these  points  [12). 

The  main  disadvantage  of  this  simple  procedure  is  very  slow  convergence  if  iteration 
is  carried  out  to  the  point  at  which  un  -  un+1  on  the  computer  [7J.  instead  of  wasting 
iterations  which  result  in  little  improvement  in  the  bounds  for  the  solution,  it  is  better 
to  stop  when  the  decrease  in  w(U"  fl  uj+1)  falls  below  a  given  limit. 

Monotonicitv  Methods 

Interval  analysis  of  integral  equations  is  simplified  considerably  if  the  function 
g(s,t,v,w)  has  known  monotonicity  properties  with  respect  to  v,w.  The  underlying  partial 
ordering  for  real  functions  is,  of  course,  componentwise.  We  say  that  v  <  w  if  v(s)  <  w(s) 
for  all  s  €  S.  If,  for  example,  g(s,t,v,w)  is  isotone  (monotone  increasing)  with  respect 
to  v,  and  antitone  (monotone  decreasing)  with  respect  to  w,  then 


G(s,t,U(s),U(t))  - 


[g(s,t,u(s),u(t) ),  g(s,t,u(s),u(t) )] 


(19) 


will  provide  an  interval  inclusion  of  g(s,t,u(s),u(t))  for  u  €  U.  Since  interval  iteration 
is  inclusion  monotone  [1],  integration  of  (19)  will  yield  an  interval  inclusion  G  of  the 
integral  operator  g.  Methods  using  monotonicity  properties  of  operators  to  obtain  lower 
and  upper  bounds  for  solutions  of  functional  equations  have  been  investigated  extensively 
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by  Collatz  [2]  and  others.  In  theory,  a  function  of  bounded  variation  can  be  expressed  as 


the  sum  of  an  isotone  and  an  antitone  function,  so  that  g(s,t,u(s ) ,u(t) )  could  be  expressed 
in  terms  of  a  finite  number  of  such  functions?  however,  it  is  not  usually  possible  to 
perform  this  construction  explicitly  [4] . 

The  actual  implementation  of  monotonicity  methods  is  hampered  by  the  fact  that  is 
difficult  to  carry  out  and  represent  the  transform  G(U)  on  the  computer,  even  for  single 
interval  functions  U(s).  For  example,  if  u  and  u  are  step  functions,  then  G(  [jj,u]  )  will 
not  in  general  be  a  step  function,  and  the  representation  of  more  and  more  complicated 
functions  will  be  necessary  as  the  iteration  proceeds.  A  solution  to  this  problem  is 
provided  by  the  concept  of  directed  rounding  in  function  space. 

Directed  Rounding 

For  actual  computation,  we  take  a  set  4  of  functions  which  can  be  represented  exactly 
on  a  computer.  For  example,  4  can  consist  of  step  functions  which  have  nodes  and  values 
which  are  floating-point  numbers.  An  upward  rounding  operator  A  from  the  set  R  of  real 
functions  into  4  has  the  following  properties: 

(i)  Au  »  u  if  u  (  4, 

(HI  u  <  Au  for  u  f  R, 

(iii)  if  u  <  v,  then  Au  <  Av. 

These  conditions  imply,  among  other  things,  that  there  is  no  element  v  of  4  such  that  u  <  v 
<  Av  [3],  hence,  this  rounding  is  "optimal"  in  a  certain  sense.  The  downward  rounding 
operator  V  is  defined  in  a  similar  way,  with  the  inequality  signs  reversed.  Thus,  if  one 
has  the  same  monotonicity  properties  for  the  function  g( s, t,u(s ) ,u(t ) )  as  in  (19),  an 
interval  inclusion  of  g  can  be  computed  by 
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G(U)  =  [V/  1G(s,t,u(s),u(t))dt, 


A/^GCSrt.uts^uCt)  )dt  ], 


(20) 


An  an  example  of  this  type  of  confutation,  consider  the  integrand 


g(s,t,u(s),u(t) ) 


lu(s)u(t) 
s  +  t 


(21) 


which  arises  in  radiative  transfer  theory,  where  0  <  X  <  0.5  [12].  Here,  t  will  be  the  set 
of  step  functions,  and  it  follows  from  (21)  that  the  integrand  is  isotone  for  positive 
functions.  Let  f  f  4  have  value  f^  in  the  subinterval  of  [0,1].  Then,  for  s  f  S^» 


g(f ) <s) 


S  +  8  . 


1  +  Xsf  •  V  f.«ln{ — —— *1— }, 
1  3  8  +  Sj-1 


(22) 


i  monotone  function  which  can  be  rounded  downward  and  upward  to  bound  the  integral  operator 
[10].  Numerical  computations  using  this  method  for  bounding  the  integral  operator  show 
noth  faster  convergence  and  applicability  to  an  extended  range  of  X  as  compared  to  the 
simple  approximation  to  the  integral  operator  using  interval  arithmetic  [10].  For  example, 
(14)  holds  when  the  operator  is  approximated  by  interval  Riemann  sums  only  for  0  <  X  < 

). 365,  while  the  use  of  (22)  extends  this  range  to  0  <  X  <  0.451  [10].  The  number  of 
Iterations  is  also  drastically  reduced  over  those  reported  in  [7]. 


Iterative  Residual  Correction 

The  method  of  interval  iteration  can  be  combined  in  a  simple  way  with  other 
tpproximate  methods  for  the  solution  of  integral  equations.  The  result  gives  essentially  a 
rathod  for  choosing  the  initial  interval  UQ  in  which  the  solution  u  of  (2)  is  sought. 
Suppose,  for  example,  that  a  numerical  solution  of  the  integral  equation  yields  the 
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approximation  uQ  »  ug(s)  to  a  solution.  If  this  solution  is  believed  to  be  correct  within 
e,  then  it  is  natural  to  start  the  interval  iteration  with  Uq(s)  =  Ug(s)  +  e*[-1/1]»  for 
example.  If  it  turns  out  that  (14)  is  satisfied  for  n  =  0,  then  the  interval  computation 
immediately  validates  the  existence  of  the  solution  u  in  the  interval  UQ,  and  also  the 
error  bound  €  13] .  The  interval  iteration  can  then  be  continued  if  it  is  desired  to  obtain 
a  better  error  bound  for  the  approximate  solution  than  e.  The  same  applies  to  Un  if  (14) 
is  satisfied  for  some  value  of  n  >  0. 

The  advantages  to  this  approach  are  that  it  may  be  much  quicker  to  obtain  a  good 
approximation  to  the  solution  of  the  integral  equation  using  floating-point  arithmetic 
instead  of  interval  arithmetic,  and  the  width  of  the  initial  interval  UQ  is  usually  very 
small.  This  is  important  because  interval  inclusions  of  functions  on  intervals  usually 
have  excess  width  over  the  range  of  the  function  on  that  interval,  that  is,  w(F(X))  can  be 
much  larger  than  w(f(X))  for  a  given  real  function  f  and  real  interval  X.  However,  it  is 
known  that  w(F(X))  ♦  w(f(X))  as  w(X)  *  0  if  f  is  continuous  [5].  This  convergence  result 
implies  that  Interval  calculations  conducted  with  intervals  of  small  width  will  estimate 
the  values  of  the  included  real  function  more  accurately,  and  the  iteration  will  hence  also 
converge  more  rapidly  to  whatever  tolerance  limit  is  set  by  the  solver  of  the  equation  or 
the  computing  equipment  used. 

In  terms  of  the  residual 


r0  “  9<u0)  "  u0' 


one  can  write 


u  "  u0  “  r0  +  g(u)  “  9(uq), 

and  if  G'  denotes  an  interval  extension  of  the  Frfichet  derivative  g*  of  g,  then  the 
interval  mean  value  theorem  [11]  gives 
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u  -  u0  €  r0  +  G’(U)(U  -  Uq), 

or,  for  O0  «  u„  +  e*[-1.1], 

u  -  uQ  €  rQ  +  e*G'(00)-t-1,1]. 

Eigenvalue  Problems 

A  simple  application  of  the  interval  method  for  solution  of  integral  equations  is  to 
the  eigenvalue  problem 

u(s)  -  X/abK(s,t)u(t)dt  (31) 

for  Fredholm  integral  operators  X  with  kernel  K(s,t).  It  follows  from  the  above  that  if 
U(s)  is  an  interval  function  such  that  0  f.  (U,U)  and  U  contains  an  eigenfunction  u  of  X, 
then  the  eigenvalue  \  belonging  to  the  eigenfunction  u  will  be  contained  in  the  interval 
Rayleigh  quotient 

A  -  (XO,U )/(U,U) ,  (32) 

where  the  scalar  products  are  formed  by  interval  integration,  and  the  division  1b  done  in 
interval  arithmetic  [8] •  This  process  forms  the  basis  of  an  interval  iteration  method. 
Starting  with  U0,  one  can  compute  Aq  by  (32).  Having  Un  and  A^,  the  next  step  is  to 
iterate 

vn+1  ”  V^n  n  vn  (29) 

to  convergence,  starting  with  V0  »  U„.  If  V*  denotes  the  converged  value  of  V,  then  one 
takes 
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and  then 


(KU  U  .) 
a  n  n±  j  n+  j 
A  1 1  .  .v  . 

"  (Vl'°n+1» 


The  interval  Rayleigh  quotient  (32)  haa  the  advantage  that  it  gives  lower  and  upper 
bounds  for  all  eigenvalues  of  K,  including  the  smallest.  Host  methods  for  approximating 
the  smallest  eigenvalue  give  only  upper  bounds  [8].  For  example,  let  K(s,t)  be  the  Green's 
function  for  the  differential  operator  -y"  on  [0,*]  with  boundary  conditions  y(0)  »  y(«)  « 


0.  One  can  start  with  UQ (s)  defined  by 


u(s)  =1,  0  <  s  <  «, 


U(s)  *  (2/ir)s,  0  <  s  <  n/2. 


U(s)  -  (2/v ) (ir  -  s),  it/2  <  a  <  *, 


and  obtain  Ag  =  [0.4052,2.5033],  while  a  second  iteration  gives  A1  -  [0.6792,1.9454],  which 
provides  better  guaranteed  lower  and  upper  bounds  for  the  smallest  eigenvalue  X  -  1  of  the 
problem. 
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